Dans le cadre de la SAE202 : Exploration algorithmique d'un problème Il nous a été demandé de réaliser un projet mathématiques et informatique. Nous avons choisis de réaliser le projet Fractales : ensemble de Mandelbrot ou Julia .
Voici la description qui a été donnée :
- Niveau de difficulté : ☆☆☆
- Description : L'ensemble de Mandelbrot est une fractale.
- Animation : Le zoom près d'un point remarquable.
En effet l'un des objectifs du projet est, entre autre, de réaliser au moins une animation. Pour ce faire nous avons entièrement réaliser le TP dédié (disponible en ligne ici).
Nous avons travaillé environ 40H sur ce projet et avons réparti le travail comme suit :
- ZHAN Pascal : Animation de l'ensemble de Mandelbrot (calcul et dessin de l'image, calcul de tableau, etc...), Programmation des fonctions pour illustrer : iterator, drawMandelBrot etc...
- JIE Juan : Recherche et définition des termes, Fonction calcul MandelBrot et getEchappement pour les couleurs.
Définition de Fractale
Une fractale est un objet mathématique, qui reste similaire à toutes les échelles.
Cet objet est "infiniment morcelé" ce qui signifie que l'on peut zoomer autant de fois que l'on veut, l'objet sera toujours entièrement représenter.
Définition d'un nombre complexe
Les nombres complexes sont une extension des nombres réels. Cette extension contient un nombre imaginaire représenter par i où i² = -1 et (-i)² = -1.
De ce fait on retrouve la formule suivante a + ib où a et b sont des nombres réels, a désigne la partie réel (abscisse) et b la partie imaginaire (ordonnée)
Cette formule permet d'écrire tous nombres complexes.
Plus précisément un nombre réel est un nombre que l'on peut placer sur un axe, l'axe des nombres réels alors qu'un nombre complexe ce place sur deux axes celui des nombres réels mais aussi celui des nombres imaginaires
Les nombres complexes restent avant tout des nombres on peut donc effectuer des opérations sur ces nombres.
Exemple d'addition
Point A = (4 + 1'i'), Point B = (2 + 4'i')
Point C = A + B = (2 + 4'i') + (4 + 1'i') = (6 + 5'i')
Multiplication des nombres réels
Pour la multiplication des nombres réels c'est un peu plus compliqué car nous avons besoin de la norme et de l'argument.
L'argument correspond à l'angle entre le nombre complexe et l'axe des réels (l'abscisse).
Pour multiplier 2 nombres complexes on commence par additionner l'argument des 2 nombres.
Suite à cela nous avons besoin de multiplier la norme des nombres complexes qui correspond à la distance entre le nombre complexe et son origine. Pour calculer cette distance nous devons effectuer la formule suivante : RACINE(a² + b²'i')
Donc, pour multiplier des nombres complexes nous avons besoin d'additionner l'argument de ces nombres et de multiplier leurs normes.
Exemple de multiplication
Point A = (2 + 3'i'), Point B = (-1 + 1'i')
--> argument = 10°
--> argument = 20°
--> Norme de A = RACINE(2² + 3²'i') = 3,6
--> Norme de B = RACINE(-1² +1²'i') = 1,4
Point C = A x B = 10° + 20° = 30°
= 1,4 x 3,6 = 5,04
Point C = (-5 - 1'i')
Maintenant que l'on a appris à manipuler les nombres complexes on va appliquer une fonction qui va nous permettre de créer un cheminement de point qui formera notre fractale, la fonction est la suivante : z --> z² - 3/4.
Point X = (1 + 1'i') --> z = argument et norme = 45° et RACINE(2)
z² = 45², RACINE(2)² = 0 + 2'i'
0 + 2'i' + 3/4 = -0,75 + 2'i'
Lorsque que l'on répète cette fonction, c'est-à -dire qu'on réapplique la formule à partir du point obtenu, on remarque que l'on obtient des nombres complexes qui s'éloigne de plus en plus de l'origine du repère.
Si l'on applique la même formule sur les nombres complexes 1 + 0,5'i' on remarque que les points s'éloigne aussi de l'origine du repère
On dira donc que les nombres 1 + 1'i' et 1 + 0,5'i' génèrent des suites non-bornées. Cependant, les nombres tels que 1 + 0,2'i' génèrent eux des suites bornées.
Voici là , les 2 catégories de suites générées par les nombres complexes : bornée et non-bornée.
Lorsque que l'on combine le calcul des suites non-bornées et le calcul des suites bornées on obtient une Fractale.
Pour être plus précis, la combinaison des suites bornées et non-bornées forme un ensemble que l'on appelle Julia et dans cet exemple on à former l'ensemble de Julia de la fonction z² - 3/4.
Nous allons avoir besoin pour la suite les bibliothèques suivantes :
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib import cm
from IPython.display import HTML
import time
L'ensemble de Mandelbrot (ou ensemble complexe de Mandelbrot) est un ensemble de points qui forment une fractale dans le plan complexe, du nom du mathématicien Benoît Mandelbrot.
Les ensembles de Mandelbrot présentent certaines similitudes avec les ensembles de Julia, comme l'utilisation du même polynôme quadratique complexe pour l'itération.
Le grossissement infini de la collection Mandelbrot peut avoir des détails subtils, et ce magnifique motif n'est généré que par une formule simple. Par conséquent, certaines personnes pensent que l'ensemble de Mandelbrot est "la figure géométrique la plus étrange et la plus magnifique jamais réalisée par des êtres humains", et qu'elle était autrefois appelée "l'empreinte digitale de Dieu".
Comme vu précedemment, le set de Mandelbrot est une collection de points qui forment un motif fractal sur le plan complexe. Il peut utiliser le polynôme quadratique complexe suivant:
\begin{equation} f_c(z) = z^2 + c \end{equation}Ici c = x + yi est un nombre complexe, si la partie réelle x et la partie imaginaire y sont considérées comme abscisse et ordonnée sur le plan complexe, alors c correspond à un point sur le plan complexe.
À partir de z = 0, le polynôme quadratique ci-dessus est appliqué à plusieurs reprises pour un calcul itératif, et la valeur de z changera continuellement, ce qui s'étendra soit à l'infini (cas a), soit restera à l'intérieur d'un cercle de rayon fini (cas b). Si c'est le cas b, on dit que le paramètre c conduit à une z-séquence itérative non divergente, auquel cas c appartient à l'ensemble de Mandelbrot, sinon ce n'est pas le cas.
L'ordinateur n'a évidemment aucun moyen d'effectuer des calculs sans fin pour voir si la valeur z itérative causée par c diverge, nous calculons donc simplement le temps d'échappement au point c - c'est-à -dire combien d'itérations se sont écoulées, la valeur z modulo supérieure à 2. Pendant le processus de dessin, ce temps d'échappement sera utilisé comme couleur du point c sur le plan.
def getEchappement(c,max_iter):
'''
c : un nombre complex
max_iter : le temps d'échappement maximal
Nous calculons le temps d'échappement pour le paramètre c qui sera utilisé comme couleur de du point
'''
z = c
for i in range(1, max_iter):
if abs(z)>2:
return i
z = z*z+c
return i
Comme indiqué dans le code, abs(z) est utilisé pour calculer le module d'un nombre complexe, où le temps d'échappement est au maximum max_iter-1 .
En même temps, il convient de noter que Python supporte nativement les opérations sur les nombres complexes, z*z + c est l'opération entre les nombres complexes, et son code n'est pas différent de l'opération entre les nombres réels ordinaires.
Ensuite, l'ensemble de Mandelbrot peut être calculé. L'ensemble de Mandelbrot généralisé est un plan complexe infini, et évidemment, nous ne pouvons en dessiner qu'un petit morceau.
Le code suivant calcule l'ensemble de Mandelbrot dans une zone carrée avec (xCenter, yCenter) comme point central et semiWidth comme demi-largeur. La plage de coordonnées x de cette zone carrée est xCenter ± semiWidth et la plage de coordonnées y est yCenter ± semiWidth.
Considérant que l'ordinateur ne peut calculer que des points discrets, la zone carrée sera divisée en N*N points pour le calcul. Lorsque N=600, la matrice résultat de la fonction suivante est une image 600x600, et la valeur de l'élément rectangulaire représente la couleur du point.
def calculMandelbrot(xCenter, yCenter, semiWidth, N, max_iter):
xFrom,xTo,yFrom,yTo = xCenter-semiWidth,xCenter+semiWidth,yCenter-semiWidth,yCenter+semiWidth
y,x = np.ogrid[yFrom:yTo:N*1j,xFrom:xTo:N*1j]
c = x + y*1j
return np.frompyfunc(getEchappement,2, 1)(c, max_iter).astype(np.float64)
y,x = np.ogrid[1:5:5j,1:3:3j]
print("y=")
print(y)
print("x=",x)
print("La forme y:",y.shape,"forme de x:",x.shape)
np.ogrid est un objet spécial qui utilise des tuples de tranche comme indices et est spécifiquement utilisé pour générer des tableaux pour la diffusion. 1:5:5j est quelque peu similaire à np.linspace(1,5,5), qui est utilisé pour générer une séquence arithmétique avec 5 éléments de 1 à 5 (y compris 5), où j ne représente pas un nombre complexe, juste une forme d'expression.
Comme vous pouvez le voir, y généré par ogrid est un tableau à deux dimensions avec 5 lignes et 1 colonne, et x est un tableau à deux dimensions avec 1 ligne et 3 colonnes. Voir ci-dessus pour sa valeur.
z = y*1j
print("z = y*1j =",z)
y*1j multiplie le nombre complexe 1j par chaque élément de y et renvoie un tableau de même forme (dimension).
Étant donné que le module numpy n'est pas le contenu natif de Python, le signe * ici n'est pas une multiplication Python normale. Il s'agit en fait d'une fonction opérateur surchargée.
Cette fonction se trouve dans le module numpy et est écrite en langage C. On peut voir que la partie réelle de l'élément résultat de y*1j est tout 0, et la partie imaginaire est égale à la valeur de l'élément d'origine multipliée par 1j.
Regardons maintenant ce que donne z + x (ici z = y*1j)
c = x + z
print(c)
print("La forme de c:",c.shape)
x est un tableau 2D avec 1 ligne et 3 colonnes, z = y*1j est un tableau 2D avec 5 lignes et 1 colonne. L'ajout des deux est en fait une fonction opérateur dans le module numpy surchargé qui ajoute les éléments correspondants des deux tableaux. La question est, deux tableaux ont des formes/dimensions différentes, comment déterminer la correspondance entre les éléments ?
Numpy effectue l'opération dite de "diffusion", et le tableau de résultat c prendra la valeur maximale de chaque axe de x et z, et la dimension résultante est (5,3), soit 5 lignes et 3 colonnes. Pour chaque élément de c, tel que c[3][2], sa valeur doit être égale à x[3][2] + z[3][2], puisque x[3][2] n'existe pas, broadcast La règle prend x[0][2] = 3 à la place, de même z[3][2] n'existe pas, et la règle broadcast prend z[3][0] = 0. + 4.j à la place, donc , c[ 3][2] = 3 + 0. + 4.j = 3. + 4.j.
Notez que la partie réelle de l'élément résultat provient du tableau x et la partie imaginaire provient du tableau z, qui est le tableau résultat de y*1j.
xFrom,xTo,yFrom,yTo = xCenter-semiWidth,xCenter+semiWidth,yCenter-semiWidth,yCenter+semiWidth y,x = np.ogrid[yFrom:yTo:N1j,xFrom:xTo:N1j] c = x + y*1jprint("c.shape:",c.shape,"x.shape:",x.shape,"y.shape:",y.shape)
Maintenant, ce programme est plus facile à comprendre. Il utilise la même méthode que le processus d'interprétation ci-dessus pour diviser la surface carrée du plan complexe avec (xCenter, yCenter) comme point central et semiWidth comme rayon en N lignes et N colonnes. Le tableau multidimensionnel c est de (N,N).
Les éléments de c sont des nombres complexes, la partie réelle varie horizontalement avec une différence égale entre xFrom et xTo, et la partie imaginaire varie verticalement avec une différence égale entre yFrom et yTo.
Regardons maintenant la dernière ligne
return np.frompyfunc(getEchappement,2, 1)(c, max_iter).astype(np.float)
getEchappement(c) est une fonction du langage Python qui a été définie auparavant. Elle accepte un paramètre complexe c et un paramètre d'itération maximale puis calcule et renvoie le temps d'échappement correspondant au paramètre complexe.
La fonction np.frompyfunc() passe les éléments complexes du tableau à deux dimensions c et le max_iter un par un à getEchappement() pour le calcul et obtient un tableau de résultat avec la même forme que c.
Combien y a-t-il d'éléments dans c, combien de fois la fonction getEchappement() sera exécutée?
Exécutez la fonction astype de ce tableau de résultats pour convertir tous les éléments du tableau de résultats en type np.float64. Les paramètres 2,1 ci-dessus indiquent que la fonction getEchappement() accepte 2 paramètres et renvoie 1 élément de résultat.
Résumé
La fonction calculMandelbrot() renvoie finalement un tableau à deux dimensions avec N lignes et N colonnes. Les éléments du tableau sont de type np.float64, indiquant le temps d'échappement du point ; le numéro de ligne et de colonne de l'élément indique indirectement que le point de l'élément est dans le résultat Coordonnées dans le plan image/complexe.
Afin d'éviter l'utilisation de variable globale, nous avons crées une classe Stockage pour stocker les variables utiles à travers les fonctions
stock ici est un objet de la classe Stockage
class Stockage:
def __init__(self):
self.c = 'flag'
self.fig = plt.figure(figsize=(8,8),dpi=100)
self.ax = plt.gca()
self.x,self.y = -0.5,0
self.sw = 1.5
self.max_iter = 700
self.pixels = 600
def drawMandelbrot(ax,xCenter,yCenter,semiWidth,N,cmap,max_iter):
"(xCenter,yCenter)-c point centre，semiWidth: le cercle，N*N pixels."
ax.set_axis_off()
#start = time.process_time()
Mandelbrot_set = calculMandelbrot(xCenter,yCenter,semiWidth,N,max_iter)
#print("Le temps de calcul est de :",time.process_time()-start,'sec')
return ax.imshow(Mandelbrot_set,cmap=cmap)
def refresh(st):
drawMandelbrot(st.ax, st.x, st.y,st.sw, N=st.pixels,cmap=st.c,max_iter=st.max_iter)
st.fig.canvas.draw()
def art(st):
return drawMandelbrot(st.ax, st.x, st.y,st.sw, N=st.pixels,cmap=st.c,max_iter=st.max_iter)
Stock_init = Stockage() --> Permet d'initialiser notre figure grâce à la classe Stockage. Cette classe contient les valeurs nécessaires pour la création de la figure.
refresh(stock_init) --> Permet de déssiner la figure en mémoire avec toutes les variables de la class Stockage.
plt.show() --> Cette ligne permet d'afficher notre figure créer avec la fonction refresh.
stock_init = Stockage()
start = time.process_time()
refresh(stock_init)
print("Le temps de calcul est de :",time.process_time()-start,'sec')
plt.show()
La fonction drawMandelbrot
Une figure dans pyplot est appelée une figure, et la figure contient un ou plusieurs axes de sous-parcelles. La fonction drawMandelbrot() est responsable du calcul et du dessin de l'image de l'ensemble de Mandelbrot sur le sous-graphe. ax.set_axis_off() masque les coordonnées du sous-image, puis calcule le tableau d'image résultant dans Mandelbrot_set de N*N pixels, puis utilise la fonction ax.imshow() pour dessiner Mandelbrot_set sur la sous-image.
Cette fonction traite les données bidimensionnelles Mandelbrot_set comme une image bidimensionnelle, et chaque élément de Mandelbrot_set correspond à un pixel de l'image. La fonction cmap dans la fonction imshow() correspond à un objet de carte de couleurs, qui est chargé de convertir la valeur de l'élément (c'est-à -dire le temps d'échappement du nombre complexe c) dans la couleur correspondante.
#Un zoom vers le point
#(0.36024044343761436323612524444954530848260780795858, -0.641313061064803174860375015179302066579494952282)
#avec un grossissement de 0,2
stock = Stockage()
stock.x, stock.y = 0.36024044343761436323612524444954530848260780795858, -0.641313061064803174860375015179302066579494952282
stock.sw = 0.2
start = time.process_time()
refresh(stock)
print("Le temps de calcul est de :",time.process_time()-start,'sec')
plt.show()
x_fin = 0.36024044343761436323612524444954530848260780795858 y_fin = -0.641313061064803174860375015179302066579494952282
x_fin et y_fin détermines les coordonnées de l'image final (La destination de notre zoom)
animMandelbrot.fig.set_figheight(7) animMandelbrot.fig.set_figwidth(7)
Cette ligne définie la taille de l'affichage de la figure.
print("\rCalculs en cours : "+str(round(num_frame/nb_frames*100, 2))+"%", " ", end='')
Cette ligne de code nous permet de savoir l'avancement de l'exécution du programme car selon le nombre de pixels, d'itération maximale et de frames l'affichage peut prendre plus ou moin de temps.
animMandelbrot.sw *= 0.7
On détermine ici la grandeur du zoom sur la figure pour chaque frame. Le zoom finale est de 0.7 puissance nb_frames
for imgNum in range(nb_frames+1):
frame = Mon_Animation(imgNum)
artist.append([frame])
La boucle a pour objectif de remplir chaque frame de notre animation.
animMandelbrot.x += num_frame/nb_frames * (abs(animMandelbrot.x-x_fin)) animMandelbrot.y += num_frame/nb_frames * (abs(animMandelbrot.y-y_fin))
Nous nous déplaçons vers le point de destination
return art(animMandelbrot)
Renvoie l'image (Ici renvoie l'image sur 1 frame)
ani = animation.ArtistAnimation(animMandelbrot.fig, artist, interval=50, blit=True)
animation.ArtistAnimation est une fonction implémenter grâce à matplotlib, cela va nous permettre d'animer la figure. La fontion à comme paramètre la figure, le nombre de frame représenter par la variable artist, l'interval en ms entre chaque frame.
HTML(ani.to_jshtml())
Cette petite ligne transforme l'affiche en HTML pour qu'on puisse agir dessus.
#Nombre de frame
nb_frames = 90
#données
x_fin = 0.36024044343761436323612524444954530848260780795858
y_fin = -0.641313061064803174860375015179302066579494952282
animMandelbrot = Stockage()
animMandelbrot.fig.set_figheight(7)
animMandelbrot.fig.set_figwidth(7)
artist = []
animMandelbrot.pixels = 1500
animMandelbrot.max_iter = 1500
#animation
def Mon_Animation(num_frame):
print("\rCalculs en cours : "+str(round(num_frame/nb_frames*100, 2))+"%", " ", end='')
animMandelbrot.sw *= 0.7
if num_frame < 8 :
animMandelbrot.x += num_frame/7 * (abs(animMandelbrot.x-x_fin))
animMandelbrot.y -= num_frame/7 * (abs(animMandelbrot.y-y_fin))
return art(animMandelbrot)
start = time.process_time()
for imgNum in range(nb_frames+1):
frame = Mon_Animation(imgNum)
artist.append([frame])
print("Le temps de calcul est de :",time.process_time()-start,'sec')
T = 19219.69
print(f'Le temps d\' exécution est de {T//60//60}h {T//60%60}min {T%60}sec')
ani = animation.ArtistAnimation(animMandelbrot.fig, artist, interval=100, blit=True)
HTML(ani.to_jshtml())